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Abstract 

Interatomic potentials of uranium dioxide are investigated on their applicability 
to model structural stabilities beyond fluorite phase by comparing with ab-initio 
results. A high pressure cotunnite phase and loosely stacking virtual crystal are 
involved in order to get a primary confidence for large-scale deformation modelings 
of UO2 under irradiation damages. The behavior of Xe atoms in UO2 fuel is studied 
with molecular dynamics simulations. Besides the ground state bubble geometry, we 
find that a planar distribution is also (meta-)stable for xenon under thermodynamic 
perturbations. The Xe atoms with a planar configuration are in a liquid state at a 
typical reactor temperature of 1000 K, which presents a modulated layer-structure 
near the interface with solid U02- A planar defects loop remains after these Xe 
atoms are released out. 
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1 Introduction 



Uranium dioxide (UO2) is a standard nuclear fuel in pressurized-water reac- 
tors. In operation in power plants or in the context of direct disposal of spent 
fuel, a clear understanding of its structural, thermodynamical and kinetical 
properties is very important. Therefore many studies have been undertaken 
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about the behavior of this material under irradiation. The crystallographic, 
elastic and thermodynamic properties, static and high frequency dielectric 
constants, point defects formation energies, activation energies for lattice mi- 
gration have been studied very well today and can be easily found in abun- 
dant literature p!|2|3|^f5|6] . A number of work on molecular dynamics (MD) 
simulations for UO2 also has been published in recent years with classical 
interatomic potentials to get various thermodynamic and transport proper- 
ties [71l8ll9lll0llll] . However, most of these studies focused on the stable fluorite 
phase, and bears a limitation that actual UO2 might deviate from this struc- 
ture locally due to irradiation damages. It is the case especially in cascade 
simulations [T2"iri3] . This provokes a question about the applicability of clas- 
sical potentials to model the behavior of UO2 with lager-scale deformations 
with defects. 

On the other hand, as a kind of major fission product, the behavior of xenon 
(Xe) atoms in uranium dioxide also has attracted considerable experimen- 
tal and theoretical attentions. Concern has been focused particularly on its 
important role in fuel swelling and this has accordingly led to a desire to ob- 
tain a greater understanding of the basic processes governing its migration 
and trapping within the fuel. The behavior of single isolated atoms of Xe 
in stoichiometric UO2 and its interaction with Schottky vacancies have been 
well studied [THI15II16|[TT] . In this paper we prefer to investigate a large-scale 
collective behavior of Xe atoms. 

From a point of view of experiment, as a kind of nuclear fuel, uranium diox- 
ide has extensive sub-grain formation (grain sub-division) under energetic fis- 
sion fragment irradiations. The restructuring starts after accumulating certain 
amount of rare gas as fission products, such as xenon in the matrix. Observing 
experimental results, one can assume that the process of sub-grain formation is 
initiated by bubble or planar structure of nanometer scale through decomposi- 
tion and reformation of U-O-Xe complex. Sub-grains develop from a repeated 
cycle of damage (ballistic excitation and local heating), recovery (relaxation) 
and quenching. To investigate this dynamical process, molecular dynamics 
(MD) study will be employed using interatomic potentials with long range 
Coulombic interactions. In this paper, we show that a planar defect on [111] 
plane in UO2 fluorite structure seems (meta-)stable when a nanometer scale 
of the defect is created. It is reasonable because the relaxation to the ground 
state bubble geometry might be prevent by an energy barrier originated from 
the rearrangement of local atomic environment. The mechanism, which we 
find by MD simulations, could be useful for understanding the nuclear fuel 
restructuring process at high burnup irradiations. 

In following sections, we compare the results of empirical shell model and pair 
potential model with first-principles calculations. The variation of relative 
energies and structural stabilities of different phases of U0 2 is studied, to 
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Table 1 

Parameters of pair potentials between Xe and UO2 





Xe-Xe (LJ) 


Xe-0 (LJ) 


Xe-U (BM) 


Xe-U (BM-fitted) 


e(eV) 


0.017 


0.0099 


2787 


4887.7 


a (A) 


4.29 


2.50 


0.41426 


0.415 



examine the transferability of these potentials beyond the standard fluorite 
phase. In section [3l we will discuss the molecular dynamics simulations of a 
planar clustering of xenon on [111] plane of fluorite UO2 in details, followed 
by a conclusion. 



2 Transferability of interatomic potential 

2. 1 Method of calculation 

Interatomic potentials are checked by calculating the energy variation of bulk 
materials with density functional theory. Plane-wave pseudopotential method 
implemented in the code VASP [r8lfT9] is used. Pseudopotentials model the 
potential yielded by nuclei and core electrons whereas the valence electron 
wave functions are expanded in a plane-wave basis. Spin-polarized local den- 
sity approximation (LSDA) with Hubbard U term correction of the energy 
functional (LSDA+U) [20.21J is employed. The parameters for Hubbard term 
are taken as U — 4.5 eV and J = 0.51 eV, which have been checked carefully by 
S. L. Dudarev et al. [2"4"f2"2T2 3j. All calculations employ projector-augmented 
wave (PAW) pseudopotentials [2"5TI2"S] with a cutoff of kinetic energy for plane 
waves of 400 eV. Integrations in reciprocal space are performed in the first 
Brillouin zone with 18 irreducible k-points for fluorite structure and at least 
28 irreducible k-points for other phases generated with the Monkhorst-Pack 
scheme [27]. The convergence of this parameter set is well checked. The energy 
tolerance for the charge self-consistency convergence is set to 1 x 10~ 5 eV for 
all calculations. 

For the part of interatomic potential models, newly parameterized shell model 
by C. Meis [7J and pair potential model by C. B. Basak [8] are used. The former 
provides the most accurate static and thermodynamic description for fluorite 
UO2 whose details are listed in the table 1 of Ref. [7J, and the latter also gives 
a reasonable result whose parameters are given in table 1 of Ref. [8]. When 
xenon is presented, a short ranged pair potential model is adopted to describe 
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interactions between Xe atom and UO2. A Lennard- Jones (LJ) model 



is employed for Xe-Xe and Xe-0 pairs whose parameters are determined by 
fitting to first principles calculated total energies of FCC-Xe and XeC>3 in 
CuAu3 phase, respectively [28]. The Xe-U pair follows Jackson's adoption p3] 
which is a Born-Meyer (BM) function 



All parameters of these potentials are listed in table [TJ 

It is necessary to point out that the Xe-U potential parameterized by Jackson 
is as a part of a shell model where he studied the Xe-UC>2 system [2]. There- 
fore the performance of this potential may have some ambiguity. We improve 
it by fitting to ab initio energies and list its parameters in tabled] as BM-fitted. 
To calculate the energy with classical potentials, GULP [29] code is employed 
where Ewald summation is involved to treat the long-ranged Coulomb inter- 
actions among uranium and oxygen atoms. To ensure the comparability of 
calculated energies, all structures have their geometry fixed both for ab initio 
and classical potential calculations. The cutoff radius of short-ranged Van der 
Waals potentials are set to 10.0 A and fixed for all calculations. 

2.2 Results and discussions 

To check the transferability of interatomic potentials beyond fluorite phase, 
cohesive energy in cotunnite phase with a space group of Pnma that stabi- 
lizes under high pressure is calculated. An experiment showed the pressure- 
induced phase transition from fluorite to cotunnite phase starts at ~42 GPa 
[5U] , whereas a first-principles calculated setout pressure is about 38 GPa [3T] . 
Local pressure due to irradiation effects is usually higher than this threshold 
value and the excess energy transferred from irradiation tracks to lattice can 
overcome the possible barrier easily. Therefore it is very possible that a portion 
of cotunnite phase of UO2 appears in irradiation damaged fuel. 

The cell volume dependence of UO2 cohesive energy is shown in figure [U 

where each unit cell contains 12 atoms (U4O8). The cohesive energy of fluorite 

° 3 

phase at a volume of 166.375 A calculated by VASP is taken as the reference 
point. The energies predicted by interatomic potential models are shifted so 
as to have the same value as ab initio results. It is evident from figure [1] 
that in high-compression and stretching regime the pair potential (Basak) 




(1) 




(2) 
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Fig. 1. (Color online) Cohesive energies of UO2 in fluorite (Fm3m) and cotunnite 
(Puma) phases calculated by ab initio (VASP) compared with those predicted 
by interatomic potential models (shell model [7] and pair potential [8]). Note the 
latter curves have been shifted to have the same reference energy as ab initio results. 



model becomes worse although a reasonable variation of energy is predicted 
around the equilibrium volume neighborhood compared with ab initio results. 
On the other hand, the shell model (Meis) provides a rather accurate energy 
variation for fluorite phase (Fm3m) within the whole studied volume region. 
However a big deviation is observed for cotunnite phase and the predicted 
energy variation is much poor by comparison with pair potential model. In this 
sense the transferability of shell model is unacceptable to treat UO2 beyond 
fluorite structure. 

A series of intermediate structures along one possible transition path from 
Fm3m to Pnma phase are also investigated"*"!. The variation of energy along 
this path is shown in figure [2j The energy of cotunnite phase has been set as 
the reference point and all curves are shifted to have the same value at this 
point. The shell model parameterized by Meis has some deviations, the largest 
deviation from first-principles result occurs at the middlemost phase where an 
energy difference of ~6.8eV for 12 atoms is obtained. The shell model has 
a large error in description structural phase transitions of UO2. Obviously 
this fault originates from the classical treatment of electric polarization in the 
semi-empirical atomic model and damaged its transferability drastically. It 
seems to be a congenital problem for shell model and we can expect a similar 
situation when apply it to other ionic materials. 



1 This path is designed according to a principle of least displacement of atoms |31| , 
which is assumed to be the most possible candidate for the true transition path. 
However, whether it is a true path or not does not affect the conclusion made here. 
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Fig. 2. Comparison of the relative energy of UO2 from Pnma to Fm3m phase 
along a designed path calculated by VASP with interatomic potential models. The 
cohesive energy of Pnma phase is set as the reference point. 

When Xe atoms congregate to a local region and form, for example, bubbles, 
the structure of UO2 near the interface boundary will departure from standard 
fluorite phase greatly and become loosely stacked. To investigate the applica- 
bility of interatomic potential models to this kind of circumstance, we employ 
a loose virtual structure in which one uranium atom occupies at (0, 0, 0) and 
two oxygen atoms occupy (1/4, 1/4, 1/4) and (3/4, 3/4, 3/4) sites in a cubic 
cell. Comparing with fluorite phase, this virtual structure can be viewed as 
derived from it by removing three uranium atoms on the face centers as well 
as the bonded oxygen atoms. Based on above discussions, we can expect a fail- 




3.0 3.5 4.0 4.5 5.0 5.5 6.0 



Lattice parameter (A) 

Fig. 3. Variation of the relative energy of UO2 in a cubic virtual structure (see text) 
along its lattice parameter predicted by VASP compared with interatomic potential 
models. The end points of these curves are set as the reference point. 



6 



ure of the shell model on this phase. The variation of cohesive energy plotted 
in figure [3] proves that it is true. The shell model parameterized by Meis [7] 
predicted a much deeper binding energy by comparing with ab initio result. 
The increment of energy under compression is slow while it is too fast under 
stretching. Oppositely, the performance of pair potential model is acceptable 
at least, although a shallower binding energy, a smaller lattice parameter and 
bulk modulus are predicted. 

When one Xe atom is inserted into the (1/2, 1/2, 1/2) site of the previous loose 
virtual cell, as figure H] showing, the performance of shell model deteriorated 
further. Its huge deviation from ab initio results bereaves its applicability to 
this circumstance. Under high compression region, the pair potential model 
also becomes poor, whose interaction of Xe-U pair is taken after Jackson that 
labeled as BM in table [TJ This potential was not well parameterized because it 
was for a shell model [H] . To improve its quality in pair potential model, we re- 
parameterize it by fitting to the ab initio cohesive energy curve. The resulted 
energy variation (associating with pair potentials of UO2 parameterized by 
Basak) is shown as the dashed line in figure HI It is necessary to point out that 
this treatment is expedient to get qualitative understanding of xenon behavior 
in U0 2 . 




4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 
Lattice parameter (A) 



Fig. 4. Variation of the relative energy of UC^Xe in a cubic virtual structure (see 
text) along its lattice parameter predicted by VASP compared with interatomic 
potential models. The end points of these curves are set as the reference point. Pair 
potentials are used to describe Xe-U02 interactions (see table [Lj) and the dashed 
curve is the result fitting to the ab initio curve to adjust the Xe-U interaction 
parameters (listed in table [T] as BM- fitted). 
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3 Classical molecular dynamics simulation 

3.1 Simulation method 



To model the collective behavior of xenon in UO2, the classical molecular 
dynamics method implemented in DL_POLY_2 [32J code is employed, as well as 
the Basak's potential set for UO2 [H] • When Xe atoms are involved, potentials 
listed in table[T]is added to the Basak's set where the parameters of (BM-fitted) 
rather than (BM) are adopted for Xe-U pairs. Previous section shows this set 
of potentials gives an acceptable description of energy variation for UC>2-Xe 
system around large scale deformations, which provides us a confidence that 
the following simulations are qualitatively correct leastwise. 

It is well known that Xe does not dissolve in UO2, which means phase sep- 
aration is inevitable in nuclear fuel when xenon is presented as a kind of 
fission product. Usually this will lead to bubbles of Xe atoms. In nuclear re- 
actor, however, Xe atoms (or clusters) may have other meta-stable structures 
around grain boundaries. Since the behavior of bubbles has been studied ex- 
tensively [T6p5] . in this paper we prefer to focus mainly on a possible planar 
distribution of xenon in UO2. 

Previous studies had suggested that Xe atoms diffuses in trivacancies (that 
is, bound Schottky trios consisting of a uranium vacancy and two oxygen 
vacancies) [33 .34.35 . 17J . On the other hand, first-principles calculation using 
pseudopotentials approach revealed that a vacancy trio with a geometry lining 
up along (111) direction is the most stable one [36J. Therefore we intermingled 
UO2 matrix by replacing bonded UO2 trios (the yellow balls in figure [5]) with 
Xe atoms on [111] planes randomly when prepare the initial configuration 
for MD simulations. [Ill] plane is chosen just because it has the largest layer 
interval distance which can mostly preserve the UO2 matrix when Xe is doped. 

It is quite clear that the conventional cell unit of fluorite phase showed in figure 
E] is inconvenient when simulating layer- related behavior along (111) direction. 
A three times larger cell unit with (111) direction aligning with Z axis is used 
instead. A 6 x 6 x 6 supercell based on this enlarged unit cell is employed 
for following MD simulations, corresponding to 648 conventional unit cells of 
fluorite structure and more than 7000 atoms. To check the size effects, we also 
use some supercells of 12 x 12 x 6 and 6 x 6 x 12, having nearly 30000 and 
14000 atoms, respectively. In all simulations, a NPT ensemble is employed 
with zero pressure and 1000 K. To remove the artificialness introduced when 
preparing initial geometries, all configurations are heated up to 3000 K and 
then cool down firstly. 5000 simulation steps are used to equilibrate system 
before gathering data. The timestep is set as 1/s; the cut-off radius is 7.5 A for 
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short-ranged Van der Waals interactions and 12 A for long-ranged Coulomb 
interaction in real-space part of Ewald summation, with a Ewald precision of 
1.0 x 10~ 5 . 



3.2 Results and discussions 



Although it is well known that bubble is the most stable geometry for Xe in 
UO2, our following simulations show that if the initial distribution of xenon is 
far from a bubble geometry and close to, for example, a planar distribution, 
then the bubble structure can not be recovered just by usual thermodynamic 
perturbations. It is perhaps due to a high energy barrier between these two 




Fig. 5. (Color online) The fluorite structure of UO2, where the yellow trios along 
(111) direction on some [111] planes are partly substituted by Xe atoms randomly 
when preparing initial structures for MD simulation. 
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Fig. 6. (Color online) Variation of number densities of Xe along (111) direction in 
UO2 simulated at 1000 K. Seven different kinds of initial distribution of Xe on [111] 
planes are considered. 
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configurations since many atoms are involved for structure rearrangement. To 
check the stability of our simulation result, the dependence on number density 
of Xe is investigated. As figure [H] shows, seven initially different configurations 
with Xe distributed on [111] planes give very similar results. Here only the 
variations of number density of xenon along (111) direction (Z axis) are plot- 
ted. The nomenclature of curves is as follows: decimal fractions denotes the 
ratio of xenon atoms to the total heavy atoms (xenon and uranium) on the 
same [111] plane, and the last letter indicates how many consecutive layers are 
involved (namely, "s" stands for single layer and "d" for double layers). For 
example [0.25_0.5_d] means two conservative [111] layers are doped, in which 
quarter uranium atoms on the first plane and half of the second plane are 
randomly substituted by Xe. Figure [6] illustrates that a layer-like structure is 
stable for a range of density of Xe in UO2. 




Fig. 7. The structure of two consecutive Xe layers in UO2 simulated at 1000 K with 
an equal amount of xenon to uranium atoms on the same [111] layer. (A) top view 
and (B) side view. 

From a point of topology, the layer structure of Xe on [111] plane of U0 2 is kept 
very well. Figure [7] gives the final equilibrium configuration of [0.5_0.5_d] after 
a 15 ps simulation in a supercell of 12 x 12 x 6 (about 8.2 x 9.5 x 5.8nm 3 ). It is 
interesting that all Xe atoms connected together to form a huge net-like cluster 
expanding across [111] plane, as the top view shows. In fact, except the case 
of [0.25_0.25_d] with the smallest number density of Xe where isolated small 
Xe clusters are obtained, all other cases produce a single large connected net- 
like cluster of Xe. The side view in figure [3(B) demonstrates a perfect planar 
distribution of Xe in UO2, and it is surprising that no any tendency towards 
bubble geometry was observed. Even increment temperature to 3500 K, quite 
similar results are still obtained. 

Since the attractive interaction between Xe atoms is very small, this kind of 
aggregation evidently drives by uranium and oxygen atoms near the interface 
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which attempt to conserve the fluorite structure of UO2 as perfect as possible. 
The Xe-UC>2 boundary actually forms a solid-liquid interface here. The radial 
distribution function (RDF) shown in figure [H] illustrates that Xe atoms are 
in a liquid state while UO2 matrix is in a solid state (see the RDF of Xe-Xe 
pairs in the inset). It can be seen clearly from figure [8] that UO2 is in a perfect 
solid state: the nearest neighbor (NN) of O-U and 0-0 pairs comprise the first 
and second peaks of the total RDF, the next nearest neighbor (NNN) 0-0 
pair and NN pair of U-U, as well as Xe-Xe NN pair contribute to the third 
peak. The total RDF of all atoms also exhibits a characteristic of solid phase. 
But the inset graph gives an essential characteristic of liquid for xenon itself. 
Therefore the behavior of Xe in UO2 is quite analogous to a system with hard 
sphere liquid in a solid. 

It is well known that near solid-liquid interfaces, for example a hard-sphere 
and hard-wall systems or a crystal and its melt, the liquid phase will present 
some oscillation of the density along the distance from the interface [37f38] . 
This behavior results from the geometrical constraining effect of the solid 
surface on the atoms or molecules of the liquid, which causes them to order 
into quasi-discrete layers, similar to the oscillation of RDF of liquid. For liquid 
Xe having an interface with FCC aluminum, similar behavior was observed 
[39] . It is very analogous to what was shown in figure [6], except here parts of 
solid fragment of UO2 are mingled with Xe liquid on the same layer. 

It would be very interesting to study with a large amount of Xe atoms 

in UO2, where some layers are occupied by xenon exclusively. We made a sand- 

7.5-r 
7.0- 
6.5- 
6.0- 
5.5- 
5.0- 
4.5- 
£ 4.0: 

CC 3.5- 
3.0- 
2.5- 
2.0- 
1.5- 
1.0- 
0.5- 

0.0 - : 

r 1 1 1 1 ' 1 ' 1 ' 1 1 r 

1 2 3 4 5 6 

Distance (A) 




Fig. 8. Radial distribution function (RDF) of U0 2 at 1000 K of [0.5_0.5_d] config- 
uration, where two consecutive [111] layers are doped with Xe by substituting half 
uranium atoms (as well as the bonded oxygen atoms to ensure the neutrality of 
simulation cell, see Figj5]). The RDF among Xe atoms plotted in the inset shows 
typical characteristics of liquid. 
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wich model that initially six consecutive [111] planes of FCC Xe are inserted 
into UO2 in a supercell of 6 x 6 x 12 (about 4.1 x4.8 x 11.7nm 3 ), to simulate the 
behavior of Xe near a solid-liquid interface. The variation of number density 
along Z axis (perpendicular to the interface) is shown in figure [HI which pre- 
dicts a layer-like (or planar) behavior of liquid Xe clearly. The layer structure 
fades away along distance from the interface, but it still maintains beyond the 
fourth layer. We also found that when UO2 cannot preserve a good crystal 
structure, the planar structure of Xe becomes blurry. Different from previous 

0.040 -r 

0.035 - 
0.030 - 

CO 

< 

£ 0.025 - 

o 
f3 

— 0.020 - 
>, 

■35 

g 0015 - 



0.005 - 
0.000- 

-60 -40 -20 20 40 

Z distance (A) 





Fig. 9. Variation of number density of Xe along (111) direction in UO2 simulated at 
1000 K. A sandwich model is employed where initially six [111] layers are occupied 
exclusively by Xe atoms, other 30 layers are perfect UO2. Inset shows the RDF of 
Xe-Xe pairs. 
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Fig. 10. Recovery of the layer structure of uranium atoms along (111) direction after 
Xe interstitials are released out from UO2 matrix, which initially has a distribution 
of xenon on [111] planes as [0.5_0.5_d], the same as figure [71 
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studies that mainly focused on geometry constrained arised from solid surface, 
our simulation suggests that the planar distribution of Xe atoms depends on 
internal state of solid phase directly (interaction among atoms) and indirectly 
(structure relaxation and surface restructure). The most interesting is that 
although different number density of xenon will leads to different interval be- 
tween layers, the interval distance of liquid Xe always equals to that of UO2. 

It is difficult to completely recover the damaged structure after Xe is released 
out at about 1000 K. This temperature is far below the melting point and 
cannot provide enough kinetic energy to overcome recrystallization barriers. 
It is quite different from the situation under irradiations, where excess energy 
from fission tracks can melt or vaporize UO2 locally and transport enough en- 
ergy to matrix to conquer the barriers. Figure [TU1 compares the layer structure 
along (111) direction of uranium before and after Xe interstitials are released 
out. In spite of in a solid state, uranium atoms within the defected region still 
rearrange the structure from 3 layers into 2 layers. However it cannot recover 
completely and many vacancies and interstitials still remained, as figure [TT1 
shown. Far from the defects left by released Xe atoms, U0 2 keeps a perfect 
fluorite structure except few oxygen interstitials. The residual defects loop 
inherits the planar geometry of foregoing Xe interstitials, partly due to the 
mechanical stability of solid UO2. The most remnant defects are uranium and 
oxygen vacancies. It is reasonable because a plenty of UO2 molecules have been 
repulsed out of the simulation cell from their positions when fission gas Xe is 
introduced. To eliminate the open space left by them, melting and annealing 
the material is necessary. 




Fig. 11. (Color online) The residual defects in U0 2 at 1000 K after Xe (see figure 
[7]) was released out: (A) is top view and (B) side view. 
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4 Conclusion 

The transferability of classical potentials of widely used shell model and pair 
potential model beyond fluorite phase was investigated by comparing their 
predicted energy variations along compression and stretching of various struc- 
tures with ab initio results. Some virtual structures and intermediate phases 
along a structural transition were involved. Results showed that despite shell 
model gives a precise description of fluorite phase of UO2, its applicability to 
other phases is questionable. The performance of pair potential model sug- 
gested that the fault might originate from the classical treatment of charge 
polarization. It arises a problem for shell model on how to control the balance 
between precision and transferability. 

Using pair potential model, we studied the planar clustering of Xe aggregating 
in UO2. Providing a quasi-layer distribution of xenon initially, our simulations 
showed it would be stable even after the melting of U0 2 . The Xe interstitials in 
U0 2 are in a liquid state and have a layer modulation of density perpendicular 
to [111] interface, which (1) has a tendency to connect all Xe atoms together, 
except the case with a low Xe density, no isolated Xe cluster was observed; 
(2) has almost the same layer distance as the solid part of UO2; (3) maintains 
beyond four layers at least. After releasing these Xe out, a planar defects loop 
left in the U0 2 matrix. We believe it has some effects on the formation of 
rim structure of nuclear fuel materials. In addition, it is necessary to point 
out that in this paper we did not consider the planar defects perpendicular to 
other direction (say, (001) or (110)) that we prefer to leave for later studies, 
where a similar conclusion can be expected. 
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